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1. Introduction 


The continental convergence (subduction/collision) normally follows the oceanic subduction 
under the convergent forces of lateral ridge push and/or oceanic slab pull (Turcotte and 
Schubert, 2002). During these scenarios, a large amount of positively buoyant materials enter 
the trench causing slow down of the convergence that, eventually, may stop. However, before 
collision ceases, convergence between the plates can continue actively for tens of millions of 
years after ocean closure as it is testified by the 50 Ma active collisions in the Western Alps 
and Himalayas (e.g. Yin, 2006). 

A remarkable event during the early continental collision is the formation and exhumation of 
high-pressure to ultra-high-pressure (HP-UHP) metamorphic rocks, which is one of the most 
provocative findings in the Earth sciences during the past three decades. Occurrences of UHP 
terranes around the world have been increasingly recognized with more than 20 UHP terranes 
documented (e.g. Liou et al., 2004), which have repeatedly invigorated the concepts of deep 
subduction (>100km) and subsequent exhumation of crustal materials (e.g. Chopin, 2003). 
It has been suggested that the HP-UHP metamorphism can be considered as a "hallmark" 
for the modern plate tectonics regime characterized by colder subduction and started from a 
Neoproterozoic time (e.g. Brown, 2006, 2007). 

The understanding of the dynamics of continental convergent margins implies several 
different but strictly correlated processes, such as continental deep subduction, HP-UHP 
metamorphism, exhumation, continental collision and mountain building. Besides the 
systematic geological/ geophysical studies of the continental convergent zones, numerical 
modeling becomes a key and efficient tool (e.g. Burov et al., 2001; Yamato et al., 2007; Gerya 
et al., 2008; Warren et al., 2008a,b; Li and Gerya, 2009; Beaumont et al., 2009; Li et al., 2011). 
The tectonic styles of continental subduction can be either one-sided (overriding plate does 
not subduct) or two-sided (both plates subduct together) (Tao and O’Connell, 1992; Pope and 
Willett, 1998; Faccenda et al., 2008; Warren et al., 2008a), as well as several other possibilities, 
e.g. thickening, slab break-off, slab drips etc (e.g. Toussaint et al., 2004a,b). Models of 
HP-UHP rocks exhumation can be summarized into the following groups: (1) syn-collisional 
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exhumation of a coherent and buoyant crustal slab, with formation of a weak zone at the 
entrance of the subduction channel (Chemenda et al., 1995, 1996; Toussaint et al., 2004b; 
Li and Gerya, 2009); (2) episodic ductile extrusion of HP-UHP rocks from the subduction 
channel to the surface or crustal depths (Beaumont et al., 2001; Warren et al., 2008a); (3) 
continuous material circulation in the rheologically weak subduction channel stabilized at the 
plate interface, with materials exhumed from different depths (Burov et al., 2001; St6ckhert 
and Gerya, 2005; Yamato et al., 2007; Warren et al., 2008a). 

In this chapter, the processes and dynamics of continental subduction/collision and HP-UHP 
rocks exhumation are investigated by the method of large-scale numerical geodynamic 
modeling. First the numerical method is described, which is followed by the numerical model 
setup and systematic thermo-mechanical numerical experiments. The discussion section 
covers a broad range of topics related to the continental subduction and exhumation. Finally 
a concluding part is presented. 


2. Numerical modeling method 


2.1 Governing equations and numerical implementation 

The momentum, continuity and heat conservation equations for a 2D creeping flow including 
thermal and chemical buoyant forces are solved: 

(i) Stokes equation 


Ox Oz = Ox (1) 


dol, Əl, oP 
= P,T 
ox j az az gp(C,M,P,T) 
where the density o depends on composition (C), melt fraction (M), pressure (P) and 
temperature (T); g is the acceleration due to gravity. 


(ii) Conservation of mass is approximated by the incompressible continuity equation 


OVyx  dvVz 
—— —— Z 2 
ox t az” A 
(iii) Heat conservation equations 
DT ) à) 
Cp (=) =- -E 4 H, +H + H, +H (3) 
oT oT 
qx = —k(C,P,T) ay 92 = —k(C, P,T) 3z 
oP fl. ys. fos 
Ha a Ta JE’ H; = Oxx Exx + 07, Ezz A- 2 On Exz 


where D/Dt is the substantive time derivative. x and z denote the horizontal and vertical 
directions, respectively. The deviatoric stress tensor is defined by Cty, o4,, lz, whilst the 
strain rate tensor is defined by éyx, xz, zz. gx and gz are heat flux components. p is the 
density. k(C, P, T) is the thermal conductivity as a function of composition (C), pressure (P) 
and temperature (T). Cy is the isobaric heat capacity. H;, Ha, Hs, Hy are radioactive, adiabatic, 
shear and latent heat production, respectively (see Table 1 for details of these parameters). 

To solve the above equations, the I2VIS code is used (Gerya and Yuen, 2003a). It is 
a two-dimensional finite difference code with marker-in-cell technique which allows for 
non-diffusive numerical simulation of multi-phase flow in a rectangular fully staggered 
Eulerian grid. I2VIS accounts for visco-plastic deformation and several geological processes 
that are described below. All abbreviations and units used in this chapter are listed in Table 1. 
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Symbol Meaning Unit 
Ap Material constant (viscous rheology) MPa "s~! 
G Cohesion (plastic rheology) MPa 
Cp Isobaric heat capacity Jk 1K! 
E Activation energy kJ mol! 
g Gravitational acceleration ms? 
G Plastic potential Pa 
Ha, Hy, Hs, Hy | Heat production (adiabatic, radioactive, viscous, latent) Wm? 
k Thermal conductivity Wm! K-1 
M Volume fraction of melt Dimensionless 
n Stress exponent Dimensionless 
P Dynamic pressure Pa 
Priuid Pore fluid pressure Pa 
Prith Lithostatic pressure Pa 
qx, qz Horizontal and vertical heat fluxes Wm-2 
QL Latent heat of melting kJ kg! 
t Time s 
T Temperature K 
Thiquidus Liquidus temperature of the crust K 
T sean Solidus temperature of the crust K 
Ve, Us Erosion and sedimentation rate ms! 
Vy, Uz Horizontal and vertical components of velocity ms 
V Activation volume J MPa`! mol`! 
X,Z Horizontal and vertical coordinates m 
a Thermal expansion coefficient x 
b Compressibility coefficient Pa~! 
Y Strain rate a 
e ij Components of the strain rate tensor s7! 
ery Second invariant of the strain rate tensor s~ 
y Viscosity Pas 
K Thermal diffusivity ms! 
A Pore fluid pressure coefficient: A = Prtuid /P Dimensionless 
u Shear modulus Pa 
p Density kgm” 
0; i Components of the viscous deviatoric stress tensor Pa 
OTT Second invariant of the stress tensor Pa 
C yield Yield stress Pa 
T Shear stress Pa 
yp Internal frictional angle Dimensionless 
X Plastic multiplier gml 


Table 1. Abbreviations and units of the variables used in this chapter. 
2.2 Boundary conditions 


For the 2D numerical models presented in this chapter, the velocity boundary conditions are 
free slip at all boundaries except the lower one, which is permeable (Burg and Gerya, 2005; 
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Li et al., 2010). Infinity-like external free slip conditions along the lower boundary imply 
free slip to be satisfied at 1000 km below the bottom of the model. As for the usual free slip 
condition, external free slip allows global conservation of mass in the computational domain 
and is implemented by using the following limitation for velocity components at the lower 
boundary: 0vx/0z = 0, 0vz/0z = —Vz/AZexternal, Where AZexternal is the vertical distance from 
the lower boundary to the external boundary where free slip (dv; /0z = 0, vz = 0) is satisfied. 
The thermal boundary conditions have a fixed value (0°C) for the upper boundary and 
zero horizontal heat flux across the vertical boundaries. For the lower thermal boundary, 
an infinity-like external constant temperature condition is imposed, which allows both 
temperatures and vertical heat fluxes to vary along the permeable box lower boundary, 
implying constant temperature condition to be satisfied at the external boundary. This 
condition is implemented by using the limitation ƏT/ðz = (Toxternal — T)/AzZexternal Where 
Texternal iS the temperature at the external boundary and AZeyternal is the vertical distance from 
the lower boundary to the external boundary (Burg and Gerya, 2005; Li et al., 2010). 


2.3 Rheological model 

A viscoplastic rheology is assigned for the model in which the rheological behaviour 
depends on the minimum differential stress attained between the ductile and brittle fields. 
Ductile viscosity dependent on strain rate, pressure and temperature is defined in terms of 
deformation invariants as: 


exp(—) 4) 


1—n 1 
n 


ductile = (en) * F (Ap) 


where ¿rp = 0.5 bij Èij is the second invariant of the strain rate tensor. Ap, E, V and n are 
experimentally determined flow law parameters (Table 2). F is a dimensionless coefficient 
depending on the type of experiments on which the flow law is based. For example: F = 
[20-#)/n] /[30+)/2n] for triaxial compression and F = 2('-2")/" for simple shear. 

The ductile rheology is combined with a brittle/plastic rheology to yield an effective 
visco-plastic rheology. For this purpose the Mohr-Coulomb yield criterion (e.g. Ranalli, 1995) 
is implemented as follows: 


Oyield = C+ Psin( Perf) (5) 
sin( pers) = sin(g) (1— A) 
— Cyield 
1] plastic = 2èm 


where yield is the yield stress. ¿ry is the second invariant of the strain rate tensor. P is the 
dynamic pressure. C is the cohesion. @¢ is the internal frictional angle. A is the pore fluid 
coefficient that controls the brittle strength of fluid-containing porous or fractured media 
(Brace and Kohlstedt, 1980). gf can be illustrated as the effective internal frictional angle 
that integrates the effects of internal frictional angle (p) and pore fluid coefficient (A). A is the 
pore fluid coefficient that controls the brittle strength of fluid-containing porous or fractured 
media. 

The effective viscosity of molten rocks (M > 0.1) was calculated using the formula (Pinkerton 
and Stevenson, 1992; Bittner and Schmeling, 1995): 


= mo expl25 + (1 — M)(~ 48 6) 
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where 1779 is an empirical parameter depending on rock composition, being 79 = 10! Pas (i.e. 
1 x 10!4 < y <2x 1015 Pas for 0.1 < M < 1) for molten mafic rocks and yo =5x 1014 Pas 
(i.e. 6 x 10! < n <8x 1016 Pas for 0.1 < M < 1) for molten felsic rocks. Successfully tested 
for a broad range of suspensions with various bubble or crystal conventions, this formula 
takes into account, other than concentration, particle shape and size distribution. 


2.4 Partial melting model 

The numerical code accounts for partial melting of the various lithologies by using 
experimentally obtained P-T dependent wet solidus and dry liquidus curves (Gerya and 
Yuen, 2003b). As a first approximation, volumetric melt fraction M is assumed to increase 
linearly with temperature accordingly to the following relations (Burg and Gerya, 2005): 


M = 0, when T < Tsolidus 


M= f= Tsolidus 


T , when Tsotidus < T < Thiquidus (7) 
liquidus ~ “solidus 


M = 1, when T > Thiguidus 


where Tsolidus and Thiguidus are the wet solidus and dry liquidus temperature of the given 
lithology, respectively (Table 3). Consequently, the effective density, peff, of partially molten 
rocks varies with the amount of melt fraction and P-T conditions according to the relations: 


Peff = Psolid — M(Psotid = P molten ) (8) 


where Psolig and Polten are the densities of the solid and molten rock, respectively, which vary 
with pressure and temperature according to the relation: 


ppt = poll — a(T — To)][1 + BCP — Po)| (9) 


where pọ is the standard density at Py) = 0.1 MPa and Tp = 298 K; « and f are the thermal 
expansion and compressibility coefficients, respectively (Tables 1 and 3). 

The effects of latent heat H; (e.g. Stitwe, 1995) are accounted for by an increased effective 
heat capacity (Cpeff) and thermal expansion («,./) of the partially molten rocks (0 < M < 1), 
calculated as 


ðM 

Cpeff = Cp + QL( 5r )P (10) 
ðM 

Aeff = a+ poe (Sy (11) 


where Cp and « are the heat capacity and the thermal expansion of the solid crust, respectively, 
and Q_ is the latent heat of melting of the crust (Table 1). 


2.5 Topographic model 

The spontaneous deformation of the upper surface of the lithosphere, i.e. topography, is 
calculated dynamically as an internal free surface by using a low viscosity (e.g. 1018 Pas), 
initially 8-12 km thick layer (thickness of this layer changes dynamically during experiments) 
above the upper crust. The composition is either "air" (1 k/m’, above water level) or "water" 
(1000 kg/m, below water level). The interface between this weak layer and the underlying 
crust is treated as an internal erosion/sedimentation surface which evolves according to the 
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Eulerian transport equation solved in Eulerian coordinates at each time step (Gerya and Yuen, 
2003b): 

a = Oz — a — Us + Ve (12) 
where Zes is the vertical position of the surface as a function of the horizontal distance vy. 
vz and vy are the vertical and horizontal components of the material velocity vector at the 
surface. vs and ve are the sedimentation and erosion rates, respectively, which correspond to 
the relation: vs; = 0, ve = Veg, when Zes < erosion level; vs = Us9, Ve = 0, when Zes > erosion 
level; where veo and vso are the imposed constant large scale erosion and sedimentation 
rates, respectively. The code allows for marker transmutation that simulates erosion (rock 
markers are transformed to weak layer markers) and sedimentation (weak layer markers are 
transformed to sediments). 


3. Numerical model design 


up to 0 km up to 4000 km 
-—— — 


Pro-continental domain —+— Oceanic domain ——_+—— Retro-continental domain 


Continental 
marginal domain 


Upper continental crust TOO"C 
Lower continental crust 
500°C 


900°C 


1300°C 


1400 1600 1800 2400 2800 3000 


1 (c) 


Fig. 1. Initial model configuration and boundary conditions. (a) Enlargement (1700 x 670 km) 
of the numerical box (4000 x 670 km). Boundary conditions are indicated in yellow. (b) The 
zoomed domain of the subduction zone. White lines are isotherms measured in °C. (c) The 
colorgrid for different rock types, with: 1-air; 2-water; 3,4-sediment; 5-upper continental 
crust; 6-lower continental crust; 7-upper oceanic crust; 8-lower oceanic crust; 9-lithospheric 
mantle; 10-athenospheric mantle; 11-weak zone mantle; 13,14-partially molten sediment 
(3,4); 15,16-partially molten continental crust (5,6); 17,18-partially molten oceanic crust (7,8). 
The partially molten crustal rocks (13, 14, 15, 16, 17, 18) are not shown in this figure, which 
will appear during the evolution of the model. In our numerical models, the medium scale 
layering usually shares the same physical properties, with different colors used only for 
visualizing plate deformation. Detailed properties of different rock types are shown in Tables 
2 and 3. 


Large scale models (4000 x 670km, Fig. 1) are designed for the study of dynamic processes 
from oceanic subduction to continental collision associated with HP-UHP rocks formation 
and exhumation. The non-uniform 699 x 134 rectangular grid is designed with a resolution 
varying from 2 x 2km in the studied collision zone to 30 x 30km far away from it. The 
lithological structure of the model is represented by a dense grid of 7 million active 
Lagrangian markers used for advecting various material properties and temperature (Gerya et 
al., 2008; Li et al., 2010). The subducting plate is pushed rightward by prescribing a constant 
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ID symbol Flow Laws E \V\n Ap ni? 
A* Air /water 0 |0/1.0/1.0 x 10-7} 1x 1018 
B* Strong wet quartzite|154] 0 [2.3] 3.2 x 1076 |1.97 x 10” 
Cc Plagioclase Anzs [238] 0 [3.2] 3.3 x 1074 [4.80 x 107 
D* Dry olivine [5328 [3.5] 2.5 x 10% [3.98 x 1016| 
E* Wet olivine |470|8 14.0] 2.0 x 10° [5.01 x 107 
F*() Molten felsic 0 |0/1.0] 2.0 x 10-9 | 5x 10! 
G*(b) Molten mafic 0 [0/1.011.0x 10-7) 1x10! 


Table 2. Viscous flow laws used in the numerical experiments. (a) yọ is the reference viscosity, 
which is calculated as 79 = (1/ Ap) x 10%”. Other parameters are illustrated in Table 1. (b) 

The molten felsic (F*) is used for partial molten sediment and continental crust. The molten 
mafic (G*) is used for partial molten oceanic crust. The rheological data in this table are from 
Ranalli (1995). 


Material) po | Cp Kb) TCI T n QL! H, [Viscousl?) |Plastic'/) 
Air 1 | 100 | 20 = : = 0 A* 0 
Water 1000/3300| 20 - - - 0 A* 0 
M,-Solid (3,4,5) [2700/1000/ k; | Ts, | Ti, [300] 2 B* 0.15 
M,-Molten (13,14,15)/2500/1500/ kı | Ts; | Tr, |300] 2 F* 0.06 
M>-Solid (6) [3000/1000/ kı | Ts; | Ti, [300/05 | C 0.15 
Mp-Molten (16) [2500/1500] kı | Ts: | Ti, [300/05 |  F* 0.06 
M3-Solid (7,8) 130001000 k> | Ts7 | Tio |380/0.25; C* 0.15 
M3-Molten (17,18) |2900/1500| ko | Tso | Tro |380/0.25|  G* 0.06 
Mg-Solid (9,10) — |3300}1000] k3 - - - |0.022 Dt 0.6 
Mg4-Molten (11) — |3200)1500] k3 - - - |0.022 E" 0.06 
Reference(s) 12) - | 3] 5 a. AAi 4 - 


Table 3. Material properties used in the numerical experiments. (a) M is used for sediment 
and continental upper crust. M? is used for continental lower crust. M3 is used for oceanic 
crust. M4 is used for lithospheric and athenospheric mantle. Numbers in the brackets are 
corresponding to material colors in Figure 1. (b) 

kı = [0.64 + 807/(Tx + 77)| exp(0.00004Pmpa), 

k2 = [1.18 + 474/(Tx + 77)] exp(0.00004Pmpa), 

k3 = [0.73 + 1293/(Tx + 77)] exp(0.00004Pmpa). (c) 

Ts, = {889 + 17900/(P + 54) + 20200/(P + 54)? at P < 1200 MPa} or {831 + 0.06P at P > 
1200 MPa}, Ts2 = {973 — 70400/(P + 354) +778 x 10°/(P + 354)? at P < 

1600 MPa} or {935 + 0.0035P + 0.0000062P? at P > 1600 MPa}. (d) T,1 = 1262 + 0.09P, 

Tr2 = 1423 + 0.105P. (e) Parameters of viscous flow laws are shown in Table 2. (f) This 
column shows the values of sin(,¢¢), which is the effective internal frictional angle 
implemented for plastic rheology. The plastic cohesion is zero in all the experiments. (Q) 
1=(Turcotte and Schubert, 1982); 2=(Bittner and Schmeling, 1995); 3=(Clauser and Huenges, 
1995); 4=(Ranalli, 1995); 5=(Schmidt and Poli, 1998). In this table, meanings of all the 
variables are shown in Table 1. Thermal expansion coefficient « = 3 x 107° K~! and 
Compressibility coefficient B = 1 x 107° MPa~! are used for all the rocks. 
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convergence velocity (Vy) in a small internal domain that remains fixed with respect to the 
Eulerian coordinate (Fig. 1). 

In the numerical models, the driving mechanism of subduction is a combination of plate push 
(prescribed rightward convergence velocity) and increasing slab pull (temperature-induced 
density contrast between the subducted lithosphere and surrounding mantle). This type of 
boundary condition is commonly used in numerical models of subduction and collision (e.g. 
Toussaint et al., 2004b; Burg and Gerya, 2005; Currie et al., 2007; Yamato et al., 2007; Warren 
et al., 2008b) and assumes that in the globally confined three-dimensional system of plates, 
local external forcing coming either from different slabs or from different sections of the same 
laterally non-uniform slab can be significant. 

Following previous numerical studies performed with similar geodynamic settings (e.g. 
Warren et al., 2008a; Li and Gerya, 2009) we design numerical models consisting of three major 
domains (from left to right, Fig. 1): (1) a pro-continental domain, (2) an oceanic domain, and 
(3) a retro-continental domain. The subducting pro-continent comprises a marginal unit and 
an interior unit. In the continental domain, the initial material field is set up by a 35 km thick 
continental crust composed of sediment (6 km thick), upper crust (14 km thick) and lower crust 
(15 km thick), overlying the lithospheric mantle (85 km thick) and subjacent mantle (540 km 
thick). The oceanic domain comprises an 8 km thick oceanic crust overlying the lithospheric 
mantle (82 km thick) and subjacent mantle (570 km thick). The material properties of all layers 
(Fig. 1) are listed in Table 3. The initial thermal structure of the lithosphere (white lines in Fig. 
1) is laterally uniform with 0°C at the surface and 1300°C at the bottom of the lithospheric 
mantle (both continental and oceanic). The initial temperature gradient in the asthenospheric 
mantle is around 0.5°C/km. 


4. Model result 


4.1 Reference model 

The reference model is designed with prescribed convergence velocity (Vy) of 5 cm/y. All the 
other configurations and parameters are shown in Figure 1 and Table 3. 

At the initial stages, the relatively strong oceanic plate subducts along the weak zone to the 
mantle (Fig. 2a). The continental margin subducts to >100 km depth, following the high-angle 
oceanic subduction channel (Fig. 2b). The significant characters are the detachment of 
subducting upper/middle crust at the entrance zone of the subduction channel with a series 
of thrust faults formed (Fig. 2b-d). A small amount of crustal rocks located in the lower 
part of the channel are detached from the plate at asthenospheric depths, indenting into the 
mantle wedge and forming a compositionally buoyant plume (Fig. 2c). Such sub-lithospheric 
plumes are discussed in detail in Currie et al. (2007) and Li and Gerya (2009). In addition, a 
partially molten plume forms in the deeply subducted oceanic plate and moves up vertically 
until it collapses at the bottom of the overriding lithospheric mantle (Fig. 2c,d). As subduction 
continues, another partially molten plume forms in the deeply subducted continental plate. It 
also moves up vertically until it collapses at the bottom of the overriding lithospheric mantle 
(Fig. 2e,f). The characteristics and 2D and 3D dynamics of this kind of plume are studied in 
detail in Gerya and Yuen (2003b) and Zhu et al. (2009). 

As continental subduction continues, partially molten rocks accumulated in the subduction 
channel extrude upward to the crustal depths (Fig. 2d,e). Then these UHP rocks exhume 
buoyantly to the surface forming a dome structure (Fig. 2f). The exhumed UHP rocks are 
mainly located near the suture zone with a fold-thrust belt formed in the foreland extending 
for about 300-400 km (Fig. 2f). P-T paths (Fig. 2, inset) show that peak P-T conditions 
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Fig. 2. Enlarged domain evolution (1300 x 600 km) of the reference model. Colors of rock 
types are as in Figure 1. Time (Myr) of shortening is given in the figures. White numbered 
lines are isotherms in °C. Small colored squares indicate positions of representative markers 
(rock units) for which P-T paths are shown (inset). Colors of these squares are used for 
discrimination of marker points plotted in P-T diagrams and do not correspond to the colors 


of rock types. 
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Fig. 3. Peak metamorphic conditions of the reference model. (a) Peak pressure condition in 


GPa; (b) Peak temperature condition in °C. 
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of the exhumed rocks are 2.5-4 GPa and 600-800 °C, respectively (also see Figure 3 for the 
peak pressure and temperature conditions of the collision zone). This indicates the UHP 
metamorphic rocks are formed and exhumed from a depth >100 km. 


4.2 Models with variable convergence velocity 

The reference model is further investigated with lower convergence velocity (2.5cm/y) and 
higher convergence velocity (10cm/y). All the other parameters are the same as in Tables 2 
and 3. 


(a) Time = 24.6 Myr 0 Er T,°c 
O 200 400 600 800 1000 


4 P, GPa > 


3 y 


4 Fates 


-g> 
1300°C 1 PA 
(b) Time = 32.6 Myr pE T, °C 
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Fig. 4. Enlarged domain evolution (800 x 225 km) of the model with lower convergence 
velocity Vy = 2.5cm/y. Colors of rock types are as in Figure 1. Time (Myr) of shortening is 
given in the figures. White numbered lines are isotherms in °C. Small colored squares 
indicate positions of representative markers (rock units) for which P-T paths are shown 
(right). Colors of these squares are used for discrimination of marker points plotted in P-T 
diagrams and do not correspond to the colors of rock types. 


In the slow convergence regime, the continental margin also subducts to >100 km depth along 
the high-angle oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 4a). 
The subducting upper/middle crust at the entrance zone of the subduction channel detaches 
with thrust faults formed (Fig. 4a-d). With continued continental subduction, partially molten 
rocks accumulated in the subduction channel extrude upward to the crustal depth (Fig. 4c,d). 
P-T paths (Fig. 4) show that peak P-T conditions of the exhumed rocks are 2-4 GPa and 
600-800 °C. 
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Fig. 5. Enlarged domain evolution (1075 x 275 km) of the model with higher convergence 
velocity Vy = 10cm/y. Colors of rock types are as in Figure 1. Time (Myr) of shortening is 
given in the figures. White numbered lines are isotherms in °C. Small colored squares 
indicate positions of representative markers (rock units) for which P-T paths are shown 
(right). Colors of these squares are used for discrimination of marker points plotted in P-T 
diagrams and do not correspond to the colors of rock types. 


In the fast convergence regime, the continental domain continues subducting along the 
high-angle oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 5a). 
Then the crustal rocks in the lower part of the channel detach from the plate at asthenospheric 
depths, intrude into the mantle wedge, and form a horizontal compositionally buoyant plume 
(Fig. 5a-c). In addition, a partially molten plume forms in the deeply subducted plate and 
moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 
5c,d), which is similar to behavior of the reference model. After the convergence ceases (1500 
km shorting, 15 Myr), the subducted continental crustal rocks in the sub-lithospheric plume 
extrude upward to the surface forming a dome structure (Fig. 5d). P-T paths (Fig. 5) indicate 
that peak P-T conditions of the exhumed rocks are 3-4 GPa and 600-800 °C, respectively. 

This parameter sensitivity studies indicate that the slower convergence produces very small 
sub-lithospheric plume (Fig. 4a), coupled subduction channel and wide collision zone (Fig. 
4d). In contrast, the faster convergence results in very large sub-lithospheric plume (Fig. 5a), 
decoupled subduction channel and narrow collision zone (Fig. 5d). Both of the models can 
obtain UHP rocks exhumation. However, the convergence velocity changes the amount of 
crustal rocks subducted to and exhumed from UHP depth (c.f. Figs 4 and 5). 
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4.3 Models with variable thermal structure of the oceanic lithosphere 
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Fig. 6. Enlarged domain evolution (900 x 225 km) of the model with higher temperature of 
the oceanic lithosphere (hot model). Colors of rock types are as in Figure 1. Time (Myr) of 
shortening is given in the figures. White numbered lines are isotherms in °C. Small colored 
squares indicate positions of representative markers (rock units) for which P-T paths are 
shown (right). Colors of these squares are used for discrimination of marker points plotted in 
P-T diagrams and do not correspond to the colors of rock types. 


The lithospheric thermal structure plays an important role on subduction/collision processes 
(e.g. Toussaint et al., 2004a, b). Therefore we investigate the sensitivity of oceanic thermal 
gradient for the reference model. In the hot model, initial thermal structure of the oceanic 
lithosphere is linearly interpolated with 0 °C at the surface (< 10km depth) and 1300 °C at 70 
km depth (compared to 1300 °C at 100 km depth in the reference model). In contrast, the initial 
thermal structure of oceanic lithosphere in the cold model is linearly interpolated with 0 °C at 
the surface (< 10km depth) and 1300 °C at 130 km depth. 

In the hot model, the continental margin subducts following the oceanic subduction channel to 
the bottom of the lithospheric mantle (Fig. 6a). Then the crustal rocks in the lower part of the 
channel detach from the plate at asthenospheric depths, intrude into the mantle wedge, and 
form a horizontal compositionally buoyant plume (Fig. 6b). The subducting upper/middle 
crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 
6b,c). With continued continental subduction, partially molten rocks in the middle channel 
extrude upward to the crustal depth (Fig. 6c,d). The subduction channel is highly coupled. 
As a result, the partially molten rocks in the sub-lithospheric plume stay at the bottom of the 
overriding lithosphere (without exhumation). 
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Fig. 7. Enlarged domain evolution (900 x 225 km) of the model with lower temperature of the 
oceanic lithosphere (cold model). Colors of rock types are as in Figure 1. Time (Myr) of 
shortening is given in the figures. White numbered lines are isotherms in °C. Small colored 
squares indicate positions of representative markers (rock units) for which P-T paths are 
shown (right). Colors of these squares are used for discrimination of marker points plotted in 
P-T diagrams and do not correspond to the colors of rock types. 


In the cold model, the continental margin subducts following the oceanic plate to the bottom 
of the lithospheric mantle (Fig. 7a). In this case, there is no sub-lithospheric plume formed. 
Consequently, the subduction channel is thicker and thicker. The subducting upper/middle 
crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 
7b,c). The subducted continental crustal rocks extrude upward to the surface forming a dome 
structure (Fig. 7c,d). The subduction channel is highly decoupled. 

The shape and characteristics of the subduction channel in the hot model (Fig. 6) is similar to 
that in the slow convergence model (Fig. 4). It indicates that both the higher temperature and 
the slower convergence can increase the rheological coupling at plate interface. As a result, 
coupled subduction channel is produced in these two models. In contrast, decoupled channels 
are formed in the colder model (Fig. 7) as well as in the faster convergence model (Fig. 5). 


5. Discussion 


5.1 Flow modes in the subduction channel 
To a first approximation, viscous channel flow can be analysed using lubrication theory (e.g. 
England and Holland, 1979; Cloos, 1982; Cloos and Shreve, 1988a, 1988b; Mancktelow, 1995; 
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Raimbourg et al., 2007; Warren et al., 2008b; Beaumont et al., 2009). Under the lubrication 
approximations, channel flow velocity is 


1 oP 
u(x) = yo eh) FU 7) (13) 


where 77 is the assumed uniform viscosity in the subduction channel, ðP /ðx is the effective 
down-channel pressure gradient, with x measured in the down-dip direction, y is the position 
in the channel measured normal to the base, h is channel thickness and U is the subduction 
velocity of the underlying lithosphere (Fig. 8a). The overlying lithosphere is assumed to be 
stationary. 


YY Mantle 
Lithosphere 


Subduction 
dominated 


From subduction 
to exhumation 


Exhumation 
dominated 


Fig. 8. Schematic diagram showing subduction/exhumation channel flow behavior in terms 
of dominating Couette (subduction) and Poiseuille (exhumation) flows (after Warren et al., 
2008b; Beaumont et al., 2009). (a) General nomenclature. (b-d) Flow types identified in the 
models. (b) Couette flow (subduction) dominates. All flow is directed downward. (c) 
Buoyant materials stagnate at bottom of channel with the Poiseuille flow effect increasing. It 
is characterized as the transition from subduction-dominated to exhumation-dominated 
channel (d) Poiseuille flow (exhumation) dominates. Buoyancy-driven exhumation starts at 
channel bottom and propagates upward. 


When nondimensional variables u’ = u/U, h! = h/H, x! = x/h and y’ = y/h are used, 
Equation 13 reduces to 
fal 412 
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where ‘ 
p — He (eP/0x) (15) 
NepfUl 
_ 2 Neff Uys 
H= (Spyax)” 116) 


where țepp, the effective viscosity and dP /ðx are averaged quantities measured at the channel- 
scale used to estimate H. This parameter H is the characteristic channel thickness for E ~ 1, 
the balance point between downward and return flows. 

Numerical models of the subduction channels can be conveniently interpreted in terms 
of the characteristic exhumation number E (Raimbourg et al., 2007; Warren et al., 2008b; 
Beaumont et al., 2009). The corresponding flow modes associated with burial and exhumation 
of UHP rocks are shown in Figure 8. The first order dynamics can be approximated 
by the above-mentioned lubrication theory for creeping flows, and characterized in terms 
of the competition between down-channel Couette flow (U(1 — y/h) in Eq. 13), caused 
by the drag of the subducting lithosphere and the opposing up-channel Poiseuille flow 
((1/2)(0P/dx)(yh — y?) in Eq. 13), driven by the buoyancy of low-density subducted 
crustal material. This competition can be expressed through the exhumation number E, 
which is a force ratio derived from the non-dimensional channel flow equation (Eq. 14). The 
actual values that determine E (Eqs 15, 16) will depend on the particular problem and its 
evolving solution. For a channel with deformable walls and no tectonic over /under-pressure, 
dP/dx = (Pm — Pc) gsin(@) where 6 is channel dip (Fig. 8a). 

Along with the characteristic E, defined at the scale of the subduction channel, space-time 
variations in the channel flow can be interpreted in terms of the local exhumation number 
E(x,t) and corresponding flow modes (Warren et al., 2008b; Beaumont et al., 2009). During 
continental subduction, E(x, t) evolves from <1 during subduction (c.f. Fig. 8b and Fig. 2a), to 
~1 during detachment and stagnation in the subduction channel (c.f. Fig. 8c and Fig. 2b,c), to 
>1 at the onset of and during exhumation (c.f. Fig. 8d and Fig. 2d-e). Buoyancy is a necessary, 
but not sufficient, condition for UHP exhumation. Among other controlling factors (Fig. 8), 
decreasing viscosity (eff) is typically most important for driving E beyond the exhumation 
threshold. In general, E(x,t) should be regarded as a measure of local exhumation potential, 
even where the local threshold value is exceeded (E>1), efficient exhumation may be impeded 
by constrictions (small 1) or high viscosities (large eff) further up the channel. 


5.2 Coupled and decoupled subduction channel 

The numerical results show that the coupled subduction channel favors lower convergence 
velocity (Fig. 4) and hotter oceanic lithosphere (Fig. 6). It is characterized by continuous 
accretion of the weak upper continental crust resulting in the development of a thick and 
broad crustal wedge. In contrast, the higher convergence velocity (Fig. 5) and colder oceanic 
lithosphere (Fig. 7) result in decoupling of the convergent plates. Transition from coupled 
to decoupled regime occurs always at the early stages of continental collision indicating that 
insertion of rheologically weak crustal material in the subduction channel is critical for the 
subsequent evolution of the collision zone (Faccenda et al., 2009). The numerical models 
confirm that HP-UHP complexes can be formed in both coupled and decoupled channels in 
the wide range of convergence scenarios (Figs 2-7). 

As discussed in Faccenda et al. (2009), coupled collision zones (which can be either retreating 
or advancing) are characterized by a thick crustal wedge and compressive stresses (i.e. 
Himalaya and Western Alps), while decoupled end-members (which are always retreating) 
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are defined by a thin crustal wedge and bi-modal distribution of stresses (i.e. compressional 
in the foreland and extensional in the inner part of the orogen, Northern Apennines). 


5.3 Thrust fault formation and exhumation of (U)HP units 


erosion 


Fig. 9. Highly-compressional regime of continental subduction (low pull force) after 
Chemenda et al. (1995, 2000). (a), (b) and (e) are successive stages of continental subduction 
in experiments without erosion. (a)-(d) show continental subduction in experiments with 
erosion. 1, overriding plate; 2, upper crust; 3, lower crust, 4, eroded material (sediments). 


One of the most important characteristics of the numerical models presented in this study is 
the formation of the thrust faults (rheological weak zones), which is followed by exhumation 
processes (e.g. Fig. 2). Similar detachment phenomenon is also documented in analogue 
models of continental subduction (Fig. 9; Chemenda et al., 1995, 1996, 2000). 

The behavior of the subducted continental crust depends on two competing effects: upward 
buoyancy and downward subduction drag as discussed in Section 5.1. Subduction drag 
within the crust and mantle drives the subduction of buoyant crustal materials into larger 
depths (Fig 2a). At the same time the buoyancy forces and also the deviatoric stresses increase 
in the subduction channel. When the materials are no longer strong enough to sustain the 
accumulated buoyancy and deviatoric stresses, the subducted continental crust will yield 
with forming the rheological weak zone (thrust fault) (Fig. 2b,c) followed by the detachment 
and exhumation of the buoyant crustal materials (Fig. 2d-f) and release of the accumulated 
buoyancy and deviatoric stresses. 


5.4 Upper crustal structure of the HP-UHP terrane 

The upper-crustal settings of many UHP terranes share a number of structural characteristics 
(Fig. 10a; Beaumont et al., 2009): (1) a dome structure cored by the UHP nappe, (2) 
domes flanked by low-grade accretionary wedge and/or upper crustal sedimentary rocks, 
(3) overlying and underlying medium- to high-pressure nappes, (4) suture zone ophiolites 
and (5) foreland-directed thrust faults and the syn-exhumation normal faults. Our numerical 
models reproduce the general characteristic upper crustal structures (Fig. 10b), especially the 
dome structure of the HP-UHP cores, the flanked and overlaid low-grade accretionary wedge, 
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Fig. 10. (a) General characteristics of UHP complexes and surrounding upper crust 
(Beaumont et al., 2009). 1, structural dome cored by UHP nappe; 2, overlying lower grade 
rocks; 3, suture zone ophiolites; 4, medium- to high-pressure nappes; 5, foreland-directed 
thrust faults; 6, syn-exhumation normal faults. Variable scale reflects size range of UHP 
complexes. (b) Structure of UHP units and surrounding upper crust in the reference 
numerical model (Fig. 2f). 


the foreland-directed thrust faults and the fold-thrust belt. The ophiolites are distributed near 
the suture zone in the reference model (Fig. 10b). 


5.5 Influence of tectonic overpressure on the P-T paths of (U)HP rocks 

The principle of lithostatic pressure is habitually used in metamorphic geology to calculate 
burial/exhumation depth from pressure given by geobarometry. However pressure deviation 
from lithostatic, i.e. tectonic overpressure/underpressure due to deviatoric stress and 
deformation, is an intrinsic property of flow and fracture in all materials, including rocks 
under geological conditions (e.g. Rutland, 1965; Brace et al., 1970; Mancktelow, 1993, 
1995, 2008; Petrini and Podladchikov, 2000). Therefore, one important question is whether 
the principle of lithostatic pressure is applicable in subduction/collision zones where 
crystallization and exhumation of HP-UHP rocks take place. Some authors have argued 
that rocks under geological conditions are too weak to support significant overpressure (e.g. 
Brace et al., 1970; Ernst, 1971; Burov et al., 2001; Renner et al., 2001; Green, 2005). Yet, Sttiwe 
and Sandiford (1994) suggested that petrologically derived P-T paths may not record depth 
changes only but stress changes. In addition, lithospheric-scale numerical models reveal 
regions where pressure may be hundreds of MPa or even several GPa higher or lower than 
lithostatic values (Mancktelow, 1993, 1995; Petrini and Podladchikov, 2000; Burg and Gerya, 
2005). The analytical solutions show that the tectonic overpressure can be as high as 60% of 
the lithostatic value in the brittle regime. In contrast, the ductile overpressure is normally 
<10% of the lithostatic value (Li et al. 2010). 

Li et al. (2010) conducted systematic numerical simulations of continental 
subduction/collision zones with variable brittle and ductile rheologies of the crust and 
mantle. In the numerical model (Figs 11, 12), the uppermost lithospheric mantle that can 
be considered as the wall of the subduction channel shows the largest tectonic overpressure 
(>1GPa and >50% of the lithostatic pressure). However, these overpressured zones rarely 
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Fig. 11. Evolution of the wedge-like subduction channel within enlarged 530 x 210 km 
domain of the original 4000 x 670 km model. Time (Myr) of shortening is given in the figures. 
White numbered lines are isotherms in °C. Small coloured squares (with ’+’ in them) show 
positions of representative markers (rock units) for which P-T paths are shown (right). 
Colours of these squares are used for discrimination of marker points plotted in the P-T 
diagrams and do not correspond to the colours of rock types. The solid curves show P-T 
paths with dynamic pressure, while the dashed curves show that with lithostatic pressure. 


or never participate in the exhumation processes. Hence, they do not influence the P-T 
conditions of geologically distributed HP-UHP rocks in nature. The main overpressure 
region that may influence the P-T paths of HP-UHP rocks is located in the bottom corner 
of the wedge-like confined channel (Fig. 11) with the characteristic magnitude of pressure 
deviation on the order of ~0.3 GPa and 10-20% from the lithostatic values (Fig. 12). The 
degree of confinement of the subduction channel is the key factor controlling this magnitude. 
The models show that the overpressure is small (~10% lithostatic) and should not affect 
in a crucial way the metamorphic mineral equilibria of the exhumed UHP rocks. The 
challenge would be to identify the geological record to actually measure precisely such 
minor deviations. An important unresolved issue concerns the question of how the 
tectonic overpressure could be potentially recorded by mineral equilibria in natural rocks. 
Thermodynamic tools used for thermobarometry of natural rocks are mainly based on 
experimental data obtained under conditions of an isotropic stress state (ie. in the absence 
of significant deviatoric stresses) and may not be directly applicable for recording dynamic 
pressure in strongly stressed rocks. Indeed, not all overpressured rocks should be strongly 
internally stressed. For example, weak (e.g. reacting, fluid rich) rock inclusions in strong 
overpressured stressed lithosphere will also have strong overpressure, but the stress state will 
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Fig. 12. The tectonic overpressure evolution of the wedge-like subduction model (Fig. 11). 
(a-d): absolute overpressure in GPa; (a’-d’): relative overpressure (ratio between tectonic 
overpressure and lithostatic pressure) in percent (%). The time (Myr) of shortening and 
tracing markers are the same as in Fig. 11. 


be isotropic, i.e. similar to the conditions of laboratory experiments. This isotropic stress 
(dynamic pressure) will be notably different from the lithostatic pressure which will be then 
directly recorded by mineral equilibria of such rock inclusions. Obviously further efforts are 
needed to experimentally study mineral equilibria in stressed rocks. 

The tectonic overpressure is also investigated for several different subduction/collision and 
exhumation scenarios (Fig. 13). In these numerical models, the overpressures are not 
significant in the mature subduction channel and/or the inner collision belt, which suggests 
that the overpressure that can possibly affect the HP-UHP rocks is mainly related to the 
wedge-like confined subduction channel (Figs 11, 12). 
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Fig. 13. The tectonic overpressure in variable subduction-collision-exhumation scenarios. 


6. Conclusions and future perspectives 


Continental collision was investigated with numerical models where an advancing 
oceanic/continental plate subducts under a fixed continent. The most important results can 
be summarized as follows: 

(1) During the processes from continental subduction to exhumation, the flow mode in 
the subduction channel changes from Couette-dominated to Poiseuille-dominated flows. 
Numerical models of the subduction channels can be conveniently interpreted in terms of 
the characteristic exhumation number. 

(2) The coupled subduction channel is formed in the models with slower convergence or 
hotter oceanic lithosphere. Whereas faster convergence or colder oceanic lithosphere result 
in decoupling of the converging plates. 

(3) The continental margin can subduct following the oceanic plate to >100 km depth, and then 
exhume to the surface forming a HP-UHP dome. The upper crustal structures of the collision 
zone in the numerical models are consistent with both the analogue model and the natural 
UHP zones. The exhumation of UHP rocks occurs in a large variety of numerical parameters. 
(4) The main tectonic overpressure region that may influence the P-T paths of HP-UHP rocks 
is located in the bottom corner of the wedge-like confined channel with the characteristic 
magnitude of pressure deviation on the order of ~ 0.3 GPa and 10-20% from the lithostatic 
values. The degree of confinement of the subduction channel is the key factor controlling this 
magnitude. The tectonic overpressure should not affect in a crucial way the metamorphic 
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mineral equilibria of the exhumed UHP rocks. The challenge would be to identify the 
geological record to actually measure precisely such minor deviations. 

The numerical modeling method is demonstrated to be a great tool to study the geodynamic 
processes in the continental convergent margins. Most of the existed numerical models in 
the relevant topics are based on the two-dimensional regimes with the along-strike variations 
ignored. However, the tectonic processes remain inherently three-dimensional. So it is quite 
significant to conduct 3D numerical modeling to investigate the dynamics of continental 
collision with applicable to the natural tectonic settings (e.g. Alpine and Himalayan collision 
belts). 
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